#!/bin/bash
#
# Copyright (c) 2014 City of Hope
# Authers: Juan Du, Dustin E Schones
# Comments and/or additions are welcome (e-mail: jdu@coh.org)
#
# Version 1.1    11/29/14

##############################################################
# ##### Please replace PATHTO with your own directory ###### #
##############################################################

Peakaboo={PATHTO}
PYTHONPATH=$PYTHONPATH:$Peakaboo/lib
export PYTHONPATH

#===============================================================
#Parameters for running Peakaboo:
# Data directory
DATADIR="."
# Data file
DATABED=""
# Data file format (can be BED3, BED6, or BEDPE)
DATAFORMAT="BED3"
# Output directory
OUTPUTDIR="Output"
# Output bed file
OUTPUTBED="Peak.bed"
# Species, for allowed species see GenomeData.py
SPECIES="hg19"
# P-value is used to define the threshold for peak calling
PVALUE="0.02"
# WINDOW_SIZE is the size of the windows to scan the genome with
WINDOW_SIZE="250"
# WINDOW_STEP is the size of steps by which the scanning window move along the genome
WINDOW_STEP="50"
# BGSIZES is the sizes of background used in the multi-scale background model
BGSIZES="10000,20000,50000"
#=================================================================

while getopts "i:f:b:o:n:s:p:w:t:g:" opt; do
    case "$opt" in
        i)  DATADIR=$OPTARG;;
        f)  DATABED=$OPTARG;;
        b)  DATAFORMAT=$OPTARG;;
        o)  OUTPUTDIR=$OPTARG;;
        n)  OUTPUTBED=$OPTARG;;
        s)  SPECIES=$OPTARG;;
        p)  PVALUE=$OPTARG;;
        w)  WINDOW_SIZE=$OPTARG;;
        t)  WINDOW_STEP=$OPTARG;;
        g)  BGSIZES=$OPTARG;;
    esac
done

if [ ! -n "$DATABED" ]; then
    echo ""
    echo "DATAFILE (-f) REQUIRED"
    echo ""
    echo 1>&2 Usage: $0 [-i DATADIR] [-f DATAFILE] [-b FILETYPE] [-o OUTDIR] [-n OUTFILE] [-s SPECIES] [-p PVALUE] [-w WINDOW_SIZE] [-t WINDOW_STEP] [-g BG_SIZES]
    echo ""
    echo "EXAMPLE:"
    echo "./Peakaboo.sh -f test-data.bed -o test-output -s mm9"
    echo ""
    exit 1
fi
    

TEMP=`expr $WINDOW_STEP - $WINDOW_SIZE`

if [ $(echo "$TEMP > 0"|bc ) -eq 1 ]; then
    echo ""
    echo "WindowStep (-t) needs to be smaller than (or equal to) WindowSize (-w)."
    echo ""
    exit 1
fi


#=================================================================
echo ""
echo "################################################"
echo "######           Peakaboo v1.1            ######"
echo "################################################"

echo "Data directory: $DATADIR"
echo "Data: $DATABED"
echo "Data format: $DATAFORMAT"
echo "Output directory: $OUTPUTDIR"
echo "Output peak file: $OUTPUTBED"
echo "Species: $SPECIES"
echo "P-value for peak identification: $PVALUE"
echo "Window size: $WINDOW_SIZE bps"
echo "Step size: $WINDOW_STEP bps"
echo "Background sizes are: $BGSIZES bps"
#=================================================================



mkdir ./$OUTPUTDIR

echo ""
echo ""
echo "Counting unique chromosomes in the data file"
cut -f 1 $DATADIR/$DATABED | sort -k 1.4,1n | uniq >> $OUTPUTDIR/chromosome-list.txt
CHROM=$(tr '\n' ' ' < $OUTPUTDIR/chromosome-list.txt)
echo "Unique chromosomes:" $CHROM

echo ""
echo "Sample from data to calculate genome-wide threshold"
python $Peakaboo/src/Get-random-coords.py -l $OUTPUTDIR/chromosome-list.txt -s $SPECIES -o $OUTPUTDIR/sample-windows.bed -w $WINDOW_SIZE >> $OUTPUTDIR/python-log.txt 2>&1

if [ $DATAFORMAT == BEDPE ]; then
  awk '{if ($9 == "-") print$1"\t"$5"\t"$3"\t"$7"\t"$8"\t"$9; else print$1"\t"$2"\t"$6"\t"$7"\t"$8"\t"$9}' $DATADIR/$DATABED | intersectBed -a $OUTPUTDIR/sample-windows.bed -b stdin -c | cut -f 4 > $OUTPUTDIR/sample-window-scores.txt
else
  intersectBed -a $OUTPUTDIR/sample-windows.bed -b $DATADIR/$DATABED -c | cut -f 4 > $OUTPUTDIR/sample-window-scores.txt
fi

THRESHOLD=$(python $Peakaboo/src/Get-threshold.py -i $OUTPUTDIR/sample-window-scores.txt -p $PVALUE)
echo "Done Calculating genome-wide threshould"
echo ""

for c in $CHROM
  do
    echo "Processing Peak Calling for" $c
    grep -w $c $DATADIR/$DATABED > $OUTPUTDIR/$c-$DATABED
    python $Peakaboo/src/Peakaboo.py -i $OUTPUTDIR/$c-$DATABED -f $DATAFORMAT -o $OUTPUTDIR/$c-$OUTPUTBED -s $SPECIES -p $PVALUE -g $THRESHOLD -w $WINDOW_SIZE -t $WINDOW_STEP -b $BGSIZES >> $OUTPUTDIR/python-log.txt 2>&1
  done

echo ""
echo "Merge the peak files"
cat $OUTPUTDIR/*-$OUTPUTBED > $OUTPUTDIR/$OUTPUTBED

echo ""
echo ""
echo "Clean up temporary files"
rm $OUTPUTDIR/chromosome-list.txt
rm $OUTPUTDIR/chr*-$OUTPUTBED
rm $OUTPUTDIR/chr*-$DATABED
rm $OUTPUTDIR/sample-windows.bed
rm $OUTPUTDIR/sample-window-scores.txt

echo ""
echo ""
echo "Done!"


